def triangular_reference_basis(x, y, basis_type, basis_index, derivative_degree_x, derivative_degree_y):
    """
    计算三角形参考单元上的基函数及其导数（基于提供的MATLAB代码实现）
    
    参数:
        x, y: 参考单元坐标（可以是标量或数组）
        basis_type: 基函数类型（201=线性基函数，202=二次元基函数）
        basis_index: 基函数索引（201类型为1-3；202类型为1-6）
        derivative_degree_x: x方向导数阶数（0,1,2）
        derivative_degree_y: y方向导数阶数（0,1,2）
    
    返回:
        r: 基函数值或其导数值（与x,y同维度）
    """
    # 处理202类型基函数（二次元）
    if basis_type == 202:
        # 零阶导数（函数值）
        if derivative_degree_x == 0 and derivative_degree_y == 0:
            if basis_index == 1:
                r = 1 - 3*x - 3*y + 2*x**2 + 2*y**2 + 4*x*y
            elif basis_index == 2:
                r = 2*x**2 - x
            elif basis_index == 3:
                r = 2*y**2 - y
            # hxm程序编号和fenicsx不一样
            elif basis_index == 6:
                r = 4*x - 4*x**2 - 4*x*y
            elif basis_index == 4:
                r = 4*x*y
            elif basis_index == 5:
                r = 4*y - 4*y**2 - 4*x*y
            else:
                raise ValueError(f"202类型基函数索引必须为1-6，当前为{basis_index}")

        # 一阶x导数（∂/∂x）
        elif derivative_degree_x == 1 and derivative_degree_y == 0:
            if basis_index == 1:
                r = -3 + 4*x + 4*y
            elif basis_index == 2:
                r = 4*x - 1
            elif basis_index == 3:
                r = np.zeros_like(x)  # 0
            elif basis_index == 6:
                r = 4 - 8*x - 4*y
            elif basis_index == 4:
                r = 4*y
            elif basis_index == 5:
                r = -4*y
            else:
                raise ValueError(f"202类型基函数索引必须为1-6，当前为{basis_index}")

        # 一阶y导数（∂/∂y）
        elif derivative_degree_x == 0 and derivative_degree_y == 1:
            if basis_index == 1:
                r = -3 + 4*y + 4*x
            elif basis_index == 2:
                r = np.zeros_like(x)  # 0
            elif basis_index == 3:
                r = 4*y - 1
            elif basis_index == 6:
                r = -4*x
            elif basis_index == 4:
                r = 4*x
            elif basis_index == 5:
                r = 4 - 8*y - 4*x
            else:
                raise ValueError(f"202类型基函数索引必须为1-6，当前为{basis_index}")

        # 二阶x导数（∂²/∂x²）
        elif derivative_degree_x == 2 and derivative_degree_y == 0:
            if basis_index == 1:
                r = 4 * np.ones_like(x)
            elif basis_index == 2:
                r = 4 * np.ones_like(x)
            elif basis_index == 3:
                r = np.zeros_like(x)
            elif basis_index == 6:
                r = -8 * np.ones_like(x)
            elif basis_index == 4:
                r = np.zeros_like(x)
            elif basis_index == 5:
                r = np.zeros_like(x)
            else:
                raise ValueError(f"202类型基函数索引必须为1-6，当前为{basis_index}")

        # 二阶y导数（∂²/∂y²）
        elif derivative_degree_x == 0 and derivative_degree_y == 2:
            if basis_index == 1:
                r = 4 * np.ones_like(x)
            elif basis_index == 2:
                r = np.zeros_like(x)
            elif basis_index == 3:
                r = 4 * np.ones_like(x)
            elif basis_index == 6:
                r = np.zeros_like(x)
            elif basis_index == 4:
                r = np.zeros_like(x)
            elif basis_index == 5:
                r = -8 * np.ones_like(x)
            else:
                raise ValueError(f"202类型基函数索引必须为1-6，当前为{basis_index}")

        # 混合二阶导数（∂²/∂x∂y）
        elif derivative_degree_x == 1 and derivative_degree_y == 1:
            if basis_index == 1:
                r = 4 * np.ones_like(x)
            elif basis_index == 2:
                r = np.zeros_like(x)
            elif basis_index == 3:
                r = np.zeros_like(x)
            elif basis_index == 6:
                r = -4 * np.ones_like(x)
            elif basis_index == 4:
                r = 4 * np.ones_like(x)
            elif basis_index == 5:
                r = -4 * np.ones_like(x)
            else:
                raise ValueError(f"202类型基函数索引必须为1-6，当前为{basis_index}")

        else:
            raise ValueError(
                f"202类型不支持的导数阶数：x={derivative_degree_x}, y={derivative_degree_y}")

    # 处理201类型基函数（线性）
    elif basis_type == 201:
        # 零阶导数（函数值）
        if derivative_degree_x == 0 and derivative_degree_y == 0:
            if basis_index == 1:
                r = 1 - x - y
            elif basis_index == 2:
                r = x
            elif basis_index == 3:
                r = y
            else:
                raise ValueError(f"201类型基函数索引必须为1-3，当前为{basis_index}")

        # 一阶x导数（∂/∂x）
        elif derivative_degree_x == 1 and derivative_degree_y == 0:
            if basis_index == 1:
                r = -np.ones_like(x)
            elif basis_index == 2:
                r = np.ones_like(x)
            elif basis_index == 3:
                r = np.zeros_like(x)
            else:
                raise ValueError(f"201类型基函数索引必须为1-3，当前为{basis_index}")

        # 一阶y导数（∂/∂y）
        elif derivative_degree_x == 0 and derivative_degree_y == 1:
            if basis_index == 1:
                r = -np.ones_like(x)
            elif basis_index == 2:
                r = np.zeros_like(x)
            elif basis_index == 3:
                r = np.ones_like(x)
            else:
                raise ValueError(f"201类型基函数索引必须为1-3，当前为{basis_index}")

        else:
            raise ValueError(
                f"201类型不支持的导数阶数：x={derivative_degree_x}, y={derivative_degree_y}")

    else:
        raise ValueError(f"不支持的基函数类型：{basis_type}（仅支持201和202）")

    return r


def triangular_local_basis(x, y, vertices, basis_type, basis_index, derivative_degree_x, derivative_degree_y):
    """
    计算物理三角形单元上的基函数及其导数
    
    参数:
    x, y: 物理坐标
    vertices: 三角形顶点坐标 (2x3数组: [x1,x2,x3; y1,y2,y3])
    basis_type: 基函数类型
    basis_index: 基函数索引
    derivative_degree_x, derivative_degree_y: x和y方向的导数阶数
    
    返回:
    基函数值或其导数值
    """
    # 提取顶点坐标
    x1, y1 = vertices[0, 0], vertices[1, 0]
    x2, y2 = vertices[0, 1], vertices[1, 1]
    x3, y3 = vertices[0, 2], vertices[1, 2]

    # 计算雅可比矩阵元素
    J_11 = x2 - x1  # x2-x1
    J_12 = x3 - x1  # x3-x1
    J_21 = y2 - y1  # y2-y1
    J_22 = y3 - y1  # y3-y1

    # 计算雅可比行列式
    J_det = J_11 * J_22 - J_12 * J_21

    # 检查单元是否退化
    if abs(J_det) < 1e-12:
        raise ValueError("退化的三角形单元，雅可比行列式接近零")

    # 计算参考坐标
    x_hat = (J_22 * (x - x1) - J_12 * (y - y1)) / J_det
    y_hat = (-J_21 * (x - x1) + J_11 * (y - y1)) / J_det

    # 根据导数阶数计算结果
    if derivative_degree_x == 0 and derivative_degree_y == 0:
        # 零阶导数（函数值）
        return triangular_reference_basis(x_hat, y_hat, basis_type, basis_index, 0, 0)

    elif derivative_degree_x == 1 and derivative_degree_y == 0:
        # 一阶x导数
        dN_dxi = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 1, 0)
        dN_deta = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 0, 1)
        return (dN_dxi * J_22 + dN_deta * (-J_21)) / J_det

    elif derivative_degree_x == 0 and derivative_degree_y == 1:
        # 一阶y导数
        dN_dxi = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 1, 0)
        dN_deta = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 0, 1)
        return (dN_dxi * (-J_12) + dN_deta * J_11) / J_det

    elif derivative_degree_x == 2 and derivative_degree_y == 0:
        # 二阶x导数
        d2N_dxi2 = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 2, 0)
        d2N_deta2 = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 0, 2)
        d2N_dxideta = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 1, 1)
        return (d2N_dxi2 * J_22**2 + d2N_deta2 * J_21**2 + d2N_dxideta * (-2 * J_21 * J_22)) / J_det**2

    elif derivative_degree_x == 0 and derivative_degree_y == 2:
        # 二阶y导数
        d2N_dxi2 = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 2, 0)
        d2N_deta2 = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 0, 2)
        d2N_dxideta = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 1, 1)
        return (d2N_dxi2 * J_12**2 + d2N_deta2 * J_11**2 + d2N_dxideta * (-2 * J_11 * J_12)) / J_det**2

    elif derivative_degree_x == 1 and derivative_degree_y == 1:
        # 混合二阶导数
        d2N_dxi2 = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 2, 0)
        d2N_deta2 = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 0, 2)
        d2N_dxideta = triangular_reference_basis(
            x_hat, y_hat, basis_type, basis_index, 1, 1)
        return (d2N_dxi2 * (-J_22 * J_12) + d2N_deta2 * (-J_21 * J_11) + d2N_dxideta * (J_21 * J_12 + J_11 * J_22)) / J_det**2

    else:
        raise ValueError("不支持的导数阶数组合")

# 以下是六个具体的基函数实现


def linear_basis_function_0(x, y, vertices):
    """线性基函数0 (1-xi-eta)"""
    return triangular_local_basis(x, y, vertices, 101, 0, 0, 0)


def linear_basis_function_1(x, y, vertices):
    """线性基函数1 (xi)"""
    return triangular_local_basis(x, y, vertices, 101, 1, 0, 0)


def linear_basis_function_2(x, y, vertices):
    """线性基函数2 (eta)"""
    return triangular_local_basis(x, y, vertices, 101, 2, 0, 0)


def quadratic_basis_function_0(x, y, vertices):
    """二次元基函数0 (1-xi-eta)(2(1-xi-eta)-1)"""
    return triangular_local_basis(x, y, vertices, 201, 0, 0, 0)


def quadratic_basis_function_1(x, y, vertices):
    """二次元基函数1 (xi)(2xi-1)"""
    return triangular_local_basis(x, y, vertices, 201, 1, 0, 0)


def quadratic_basis_function_2(x, y, vertices):
    """二次元基函数2 (eta)(2eta-1)"""
    return triangular_local_basis(x, y, vertices, 201, 2, 0, 0)
